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ABSTRACT 


A Monte Carlo computer model was developed to simulate the 
propagatıon of light through a scattering/absorbing medium 
using various parameters and phase functions. The model 
permits characterization of the spatial and temporal spread of 
light traversing plane-parallel clouds. t was found that 
both the time and spatial spread of light in a scattering 
medium are independent of the details of the phase function for 


a cloud thickness of greater than 15 extinction lengths. 
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L. INTRODUCTION 


As HISTORY 

The existence of clouds and fog in many regions of the 
earth presents a formidable problem to the designer of an 
optical communication system whose transmission channel is 
the atmosphere. A scattering medium can inhibit system per- 
formance by inducing beam spread, dispersion in angle-of- 
arrival, degradation of spatial coherence and dispersion in 
time and frequency of the modulated optical beam. The develop- 
ment of an optical communication system for atmospheric 
applications requires an accurate knowledge of the effects of 
scattering on light propagation. Numerous studies have been 
made in an attempt to provide this knowledge. Different 
aspects of the problem have been assaulted using various 
mathematical models. Monte Carlo computer simulation has 
been used by Bucher [1], Plass and Kattawar [2] through [5], 
Junge [6] and Danielson, Moore and van de Hulst [7], to mention 
only the few used extensively in this work.  Equally prominent 
attempts have been made using analytical methods such as those 
by Arnush [8], Gordon [9], Stotts [10], Hansen [11], Ishimaru 
[12], Kennedy [£13] and Lutomirski and Yura [14]. Complicated 
numerical techniques have been used only by Dell-Imagine [15] 
and Zachor [18]. Each attempt has contributed to the stock- 
pile of knowledge required but additional problems remain to 


be investigated. Many books have been published on the subject 
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of light scattering. Those found most useful to this work 
were van de Hulst [18], Kerker [19], Chandrasekhar [20] and 
Deirmendjian [21]. Interested parties and reasons for their 


interest are briefly discussed in the following paragraphs. 


Be PURPOSE 

Current navy operational communications systems suffer 
from a umbes ог problems. There is no operational system 
which is not subject to jamming, intercept, spoofing and 
direction finding. In addition, it appears that optical 
communications systems have great promise in solving these 
problems for many applications [17]. Any information that 
could be used in evaluating such systems is desired. 

The Navy's effort to create a worldwide satellite-to-sub- 
marine communication network is another rapidly progressing 
area requiring information of the nature being investigated 
in this thesis. In this situation, light propagation through 
water further complicates the matter. 

Among other parties that may require information on light 
propagation are the United States Coast Guard in their aids 
to navigation system, NASA in satellite communications and 
researchers trying to determine the composition of the atmos- 
Phere by analyzing scattered radiation. 

The purpose of this work is to characterize both the 
Spatial and temporal aspects of a light beam propagating 
through a scattering/absorbing medium using an analytically 


verified Monte Carlo model of the system. Models were created 
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and tested against each other until an adequate level of 
confidence in each was obtained. The models were then used 
to study more specific questions concerning time and spatial 
spread. More explicitly, the models were used to observe 
what effect different phase functions had on .the temporal and 
spatial aspects of light scattering. It was believed that 
the pronounced forward peak of some phase functions would 
generate different time and spatial character than a moder- 
ately peaked phase function. An investigation attempting to 
verify this belief was an important part of this work. 

Spatial characterization of light passing through a cloud 
of finite thickness, and transition from the forward scatter 
to diffusion region was also investigated. 

Once confidence was established in the models used, the 
general areas of agreement or disagreement were studied in 
an attempt to specify when one model may be used more effici- 
ently than another for generating information concerning a 
Specific situation. The drawbacks and strong points of each 


model type have been determined and mentioned where appropriate. 
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Il. THEORETICAL PRELIMINARIES 


A. INTRODUCTION 

The purpose of this section is to review the various 
theoretical formulations currently used to characterize 
optical propagation in a single and multiple scattering medium. 
A complete coverage of any one of the topics in the following 
sub-sections is most certainly not contained here but suffici- 
ent references are provided for the reader inclined to 
investigate any topic in more detail. Only those details 
related closely to this work are described in any detail and 
even then only the important conclusions are mentioned in 
many cases. 

The order of presentation is as follows. The principal 
characteristics of the scattering/absorbing medium are men- 
tioned briefly, thereby defining frequently used terms. A 
brief summary of Mie theory is included to present the valuable 
insight necessary in investigating the problems posed. 
Radiative transfer theory is mentioned in passing because the 
solution of the complicated radiative transfer equation is 
essentially the subject of many analytical attempts at solving 
the scattering problem. This section then probes into mathe- 
matically modeling the geometry of a real atmosphere. Three 
general categories of models are discussed. 

Analytical models, which made various assumptions allowing 


closed form expressions for certain characteristics of the 
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medium, are discussed. The approximations in such models 
normally detract from the total worthiness of such an 
approach. 

The discussion then progresses into the use of numerical 
methods in evaluating analytical solutions as well as in 
evaluating the radiative transfer equation directly. This 
often overlooked powerful technique has been used by theore- 
ticians because of its ease in solving complicated integro- 
differential equations for which exact solutions have been 
proven impossible. 

Monte Carlo methods, used predominantly in this work, are 
then discussed. This technique is used extensively due to its 
ease in adaptation toodd geometries of the scattering medium. 
In this work this benefit was used only to model finite 
clouds in an otherwise homogeneous atmosphere. 

This section is completed with a general discussion of 
the results predicted by theoretical insight into the problem 


of light transfer in a random medium. 


B. PRINCIPAL CHARACTERISTICS OF SCATTERING/ABSORBING MEDIA 
This section divides the characteristics of a scattering 
medium into three sub-sections. The first defines the para- 
meters of the atmosphere which during the course of this work 
were considered constants. The second considers the types of 
scattering encountered in a scattering medium and the third 
briefly describes the significance of the phase function in 


scattering theory. 
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l. Parameters 
In general, a medium may exhibit both scattering and 
absorption. Some of the energy of a collimated incident beam 
on a particle will be scattered over all possible directions; 
the rest will be absorbed by the particle and lost from the 
radiation field. In the single scattering problem the inten- 
sity of radiation at any point along the direction of plane wave 


propagation is given by: 


-(K 4 IR -(K , *K IR 
I(R) = I(O)e S = 1(0)е E o MEN dee 


` 


where the absorption loss and scattering loss together is 
called extinction. The ratio of scattering coefficient to 
the extinction coefficient is called the albedo of single 


scatter. Thus, 


K 
de = gp AH << © « 1 (2) 
abs sca 


It is sometimes more convenient to discuss optical dimensions 
in terms of mean free paths of photons, i.e., the average 
distance between collisions. The mean free path is often 
called an extinction length and is given by the reciprocal of 
the extinction coefficient. Optical thickness of a given 
medium is a dimensionless number representing its thickness 
in extinction lengths or mean free paths. 

The transmittance of the atmosphere is the fractional 


intensity remaining in a beam after traversing a path R units 
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long. Explicitly, with K,,, in km"! 


Ме? апа К іп kılometers, 


) - e š (3) 


The term visibility is defined as [22] 





ү = |; Т 1 = 33912 I (4) 
K .02 K 
ext ext 





thus the transmittance for a path length just equal to the 

visibility is two per cent. In many atmospheric light propaga- 

tion problems the albedo is close to unity and Km in 

the above equations can be expressed by КЕСЕ with little error. 
In a real atmosphere there is a contribution to scat- 

tering by small particles as well as a contribution by large 

particles. ЕСЕ then is defined in even more detail. The 


following paragraphs present some of these details. 
2. Types of Scattering 
Atmospheric light scattering can be classified in two 
general categories. Small particle scattering, where the 
particle radius is much scaller than the wavelength of the 
incident light, is called Rayleigh scattering, and other 
scattering from larger particles is termed Mie scattering. 


The terminology of this work refers to Kea as the portion of 


y 
the scattering coefficient due to Rayleigh scattering and 


KMie as the portion due to Mie scattering. Therefore, 
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K = K + Ku T K (5) 


Mie scattering is often called particulate scattering while 
Rayleigh is often referred to as molecular scattering. A 
ratio of particulate to total scattering is defined for use 


in this thesis as: 


O i 


Scattering may be isotropic (scatters in all direc- 
tions equally) or anisotropic (scatters as a function of 
angle off incident direction), The latter only is used here. 
The problem of defining scatter for single particles is well 
understood. The theory of scattering has been extended to a 
medium of many equal sized particles, a monodispersion, and 
to a medium of many different sized particles, a polydisper- 
Exon £18, 21]. 

When a photon undergoes only one collision in tra- 
versing a scattering medium, single scattering has taken place. 
On the other hand, when a photon undergoes more than one 
collision, multiple scattering has occurred. This work in- 
volves independent single and multiple scattering in spherical 
polydispersions, where no scattering event affects other 
scattering events. For further details related to scattering 


parameters, types of scattering and other more specific topics 
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the reader is referred to an excellent text on the subject, 
McCartney [22]. The following paragraphs focus on the angular 
dependence of scattering for scattering of different types. 
3. Phase Function 

The phase function expresses in a formal manner the 
angular dependence of scattering. The phase function, denoted 
here by P(9), is defined by van de Hulst [18] as the ratio of 
the energy scattered per unit solid angle in a given direction 
to the average energy scattered per unit solid angle in all 
directions. This definition requires that the integral of 


the phase function be normalized to unity, which is to say that, 


T 
2 P(8) sin 9 dedo = 1 (7) 


When expressed in this way the phase function is a scalar 
function. In a more complex analysis of scattering theory, 
there are phase functions of different types for different 
polarizations of light. These functions form elements of the 
normalized scattering matrix which is used to transform 
Stoke's parameters. Detailed discussion of the parameters is 
given by Refs. 18 and 21. In Appendix B a program to evaluate 
the Stoke's parameters for the spherical polydispersions used 
in this work is presented and explained. The scalar phase 
function is obtained by averaging the two Stoke's parameters 
obtained using this computer adapted Mie theory. 

In many cases the phase function is easily approximated 


as a closed form function of scatter angle. In these cases, 
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the Henyey-Greenstein function is often used as the approxi- 
mating function because of its simple form. For the Henyey- 
Greenstein phase function (Figure 1), 


2 
P(8) = — = , -1«G«1 (8) 


47 (1*6?-28c0s0) 1? 

where G is a parameter which equals «cos0» for the phase func- 
tion. A G value between .80 and .87 approximates real 
atmospheric phase functions quite well in most cases. However, 
the following phase functions were not represented well using 
the ey Eresnstedn phase function no matter what parameter 
was used. 

Phase functions used in this work were of Water Cloud 
C.2 [21] at a wavelength of .53 microns (Eier 2) which was 
calculated using the program of Appendix B, and NOSC Fog 
(Figure 3) which was calculated in a similar manner using an 
experimentally determined ТК distribution near San Diego, 
California. For further discussion of phase functions, see 
"Methods of Simulation" later in this thesis. References 18 
and 21 study the phase function and its generation extensively. 

In general the phase function is quite symmetric in 
the forward and backward direction for particles small compared 
to wavelength. For large particles the phase function has an 
increasingly complicated dependence on angle of scatter. 

As one might imagine, the scalar phase function plays 


a significant role in single and multiple scattering problens. 
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Much attention has been focussed on defining the conditions 
under which the entire Mie series can be approximated by 
analytic expressions such as the Henyey-Greenstein function 
mentioned previously [l1, 7, Il, 18]. This thesis is part of 


the quest to define these conditions. 


C. MIE SCATTERING 

Mie theory considers the scattering and absorption which 
occurs when an electromagnetic wave is incident on a spherical 
particle. It generates desired parameters of the atmosphere, 
namely the extinction and scattering cross sections, as well 
as the phase function. Mie theory is extremely versatile in 
that it can be applied at any particle size to wavelength 
ratio and can be extended to many particles and more important- 
ly polydispersions. Appendix B explains the adaptation of Mie 
theory to machine computation for collecting scattering data 
for polydispersions of a given particle distribution. As 
mentioned before, Mie theory was used to generate optical 
parameters of the Water Cloud C.2 [21] model at a wavelength 
of .53 microns. It is prudent at this time to define the Mie 


Scattering parameter as 


x s E ES) 


where r is the radius of the particle and A is the wavelength. 
For the reader interested in complete derivations of 
Mie theory, see Chandrasekhar [20], Sekera [23], Born and Wolf 


[24] or van de Hulst [18]. 
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D. RADIATION TRANSFER 

Any discussion of light propagation through the atmosphere 
would not be complete without at least a mention of the 
radiative transfer equation. This equation models mathemati- 
cally what 1s happening as light transverses a scattering 
medium. The radiative transfer equation in its simplest form 
ís [10], 


d = = 
(u e POIL Tcp ну бу Т) = QN DUE CINQ E) 


(10) 

27 1 = 
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where I is the scattered radiance, I. is the radiance due to 
the distributed source produced by the incident beam transvers- 


ing the region, 4 = cos®, т = K, z ls the optical thickness, 


xt 
» are the spatial coordinates, 0, 0' are angular coordinates, 
Po is the albedo, P is the scalar phase function and t is time. 
Approximations can be made to solve this equation for 

T(t, Y, u, 8, t) in closed form. The equation has basically 
two components. One represents the field caused by the 
incident spatially distributed field and the other represents 
the field scattered out of the direct beam but redirected back 
into the same direction by other scattering events. For a 
complete derivation and explanation of this equation, 
Chandrasekhar [20] is suggested. Under specific conditions, 
approximate solutions to this equation can be found. The 


next section mentions a few of the efforts made toward finding 


useful solutions. 


Ze 





E. ANALYTICAL METHODS 

As mentioned previously, several authors have attempted 
to provide knowledge of the effects of particulate multiple 
scattering on light propagation by applying analytical methods 
[8-13, 20, 25, 26]. One analytical development by Arnush [8] 
utilized radiative transfer theory and the small angle 
approximation to characterize the light.  Arnush made two 
assumptions that greatly weakened his model: (1) He assumed 
the incident signal did not experience pulse broadening while 
traversing an optically thick medium, and (2) he assumed 
the incident beam never directly created internal emission 
sources. Stotts [10] does consider the previous details but 
still makes tae ema angle approximation,claiming that the 
phase function is highly peaked in the forward direction. 
Ine РАА at gathering knowledge include that by Ishimura 
and Hong [12] who reported an analytical study of coherence 
using first order approximations. The results are valid for 
weak fluctuations in the mediun. 

Zachor [16] uses a double-integral transform method which 
is evaluated recursively to obtain the aurole radiances contri- 
buted by successive scattering orders. He has assumed in his 
calculations a homogeneous unbounded atmosphere and his results 
are good for short ranges only. 

In any case, because of the complexity of the integrals 
involved , exact analytical methods yield results only after 


repeated numerical integration. Even in the small angle 
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approximation, the closed form solution [27] requires seven 
successive integrations to obtain a single value of the 
irradiance caused by a unidirectional point source of light. 
All this numerical work, besides being tedious, tends to mask 
the functional dependence of the results on the underlying 
physical geometry and optical parameters. 

There are many complications in solving the multiple 
Scattering problem analytically. Nevertheless, analytical 
work does contribute significantly to overall knowledge of 
the problem because in many cases the approximations made do 
not substantially deviate from reality. 

In this thesis, two of the more straightforward analyti- 
cal methods are used; Gordon [9] and Stotts [26]. The former 
concerns spatial spread and the latter, time spread. 

Dell-Imagine [15] derives a solution to the radiative 
transfer equation analytically in the usual fashion using no 
assumptions until he has to actually get numerical results. 
His numerical approach to evaluating the transfer equation 
Seems quite attractive and is briefly mentioned in the follow- 


ing section. 


F. NUMERICAL METHODS 

The solution of the radiative transfer equation is suffi- 
cient to specify the properties of a received signal which has 
passed through a multiple scattering region. The solution, 
however, requires numerical computation on a digital computer. 


The solution has the general form [15], 


ou 





where, 
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0. = scattering efficiency factor с = speed of light 


which has not been restricted to any shape of cloud. By 


establishing boundary conditions, initial conditions and 


density of particles n, a solution for I(x,y,z,t,0,0) can be 


obtained by approximating the equation by a group of simul- 


taneous algebraic equations.  Dell-Imagine describes this 


approximation and the numerical methods used in solving the 


equation and the reader is referred to his work for further 


details. 


G. MONTE CARLO METHODS 


The Monte Carlo method can be applied to any problem if 


one knows the probability for each step in the sequence of 


3l 





events and desires the probability of the total of all possi- 
ble events. Thus, the Monte Carlo method may be used to 
study problems in radiative transfer. As mentioned in the 
introduction, many authors have described their attempts to 
do so. The Monte Carlo calculations are relatively easier to 
model than other methods especially when the geometry becomes 
difficult. The Monte Carlo method is very consumptive of 
computer time and only approximate information is obtained. 
Some generally useful results have been presented by Bucher 
[1], Junge [6] and Hansen [11]. Monte Carlo results generally 
yield excellent agreement with experimental data, however, 
the results deteriorate statistically at large ranges or 
narrow observation angles. 

Details of a program for generation and curve fitting of 
Monte Carlo data are contained in Junge [6]. This thesis 
uses the groundwork of that report to extend study into 


regions other than ultraviolet light in a homogeneous space. 


=Æ THEORETICAL RESULTS 

The purpose of this section is to summarize the effects 
scattering has on light propagation through clouds as pre- 
dicted by the theories of the previous sections. In doing 
this the definitions of effects investigated are given. 
Actual quantitative relations used in this thesis are reserved 
for later sections and qualitative descriptions are found 


here. 
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Previous Monte Carlo, analytic, and numeric studies of 
the multiple scattering problem have found that light pulses 
are distorted in the following manner: 

l. Angular Spreading 

Angular spreading constitutes an effective decollima- 
tion of the incident radiation. It contributes to beam 
Spread, dispersion in angle-of-arrival and loss of spatial 
coherence. 

2. Spatial Spreading 

Spatial spreading indicates the dimensional increase 
of the beam's finite cross section. It is related to the 
above parameter in that angular spread inside the medium 
produces spatial spread as the pulse propagates on. Spatial 
coherence and beam spread are related to this parameter. 

3. Multipath Time Spreading 

Different distances along various possible propagation 
paths imply different transit times for a photon. Thus, a 
short optical pulse will incur pulse broadening after tra- 
versing a multiple scattering region. This pulse broadening 
is called multipath time spreading. The maximum pulse frequ- 
ency that can be used in communicating is determined by this 
Parameter. 

Multipath time spread is often defined as the differ- 
ence between the average transit time incurred from multiple 
scattering and the normal transit time in the absence of 
scattering. This definition is normally used when the time 
dependence of the input is that of a delta function. The 
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previous description of multipath time spread concerns input 
pulses of finite width. 
4. Total Transmission 

Total transmission describes the amount of irradiance 
left in the pulse after traversing the medium. It is in- 
versely related to the beam attenuation. 

Bucher [1] found that the amount and distribution 
of multipath time spreading was essentially independent of 
the detailed shape of the scattering function for sufficiently 
thick clouds. He also observed that the propagation para- 
meters for sufficiently thin clouds were dependent both on the 
cloud parameters and on the scattering function. 

Gordon [9] found that within certain ranges and angles 
of practical interest, the flux and beam spread function can 
be adequately approximated by closed form expressions. His 
work was related to scattering under water but can be applied 
equally well to atmospheric scattering when the phase function 
is very forward peaked. 

Time spreads on the order of microseconds to milli- 
seconds have been reported by Bucher [1], Stotts [10,26] and 
Ishimura and Hong [12]. Danielson, Moore and van de Hulst [7] 
and Hansen [11] found that the Henyey-Greenstein function 
adequately approximated the true cloud or haze phase function 
in determining the reflection and transmission characteristics 
of clouds. 

It is the purpose of this thesis to further verify 


results found previously and to quantify some of the more 
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important cloud characteristics with respect to time and 


spatial spread. 
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III. STATEMENT OF’TEE PROBLEMS 


A. INTRODUCTION 

There are many areas in the study of multiply scattered 
light that could easily be subject to further study. It was 
necessary to select a few problems which were assailable using 


available groundwork and techniques. 


EI THE PROBLEMS 

The Monte Carlo routine of Ref. 6 was employed to investi- 
gate the following problems: 

l. Dependence of Spatial Spread on Phase Function 

Characterize the spatial spread of light traversing 

Scattering media of various types. Select phase functions of 
high, moderate and low forward peakedness all with the same 
(cos0» and the same root mean squared scatter angle. Deter- 
mine the effect on spatial character of light of using the 
different phase functions. Determine the conditions concern- 
ing dependence on phase function. 

2. Dependence of Time Spread on Phase Function 
Characterize the time spread of light traversing 
Scattering media of various types. Use the same phase functions 
Selected above to determine the effect on temporal character of 
light of using the different phase functions. Determine the 


conditions concerning dependence on phase function. 
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3. Transition Region from Forward to Multiple Scattering 


Using the results of the previous two characterizations, 
determine information concerning the transition from the region 
of primarily forward scatter to the diffusion region where 
scatter is directed in all directions. Estimate the accuracy 


of the information and verify its agreement with other theories. 


4. Spatial Characterization of Light Transversing Finite 
Clouds 


Characterize the spatial spread of light traversing 
a homogeneous medium.  Characterize the spatial spread of light 
traversing the same medium except that a dense cloud of 
Scatterers of finite thickness has been added. Compare the 
two characterizations and discuss the results. 

5. Model the Above Problems 

Create a model to simulate the geometries of the above 

problems and verify at each step that the model is indeed 


generating results consistent with existing theory. 


37 





IV. METHODS OF SIMULATION 


A. MONTE CARLO METHOD 

The computer routine of Ref. 6 was avallable for use at 
the outset of this endeavor. The routine could calculate 
both spatial and temporal characteristics of light travers- 
ing a scattering medium. It required as input the scattering 
parameters and phase function parameters of the scattering 
particles of the medium. The routine was restricted to homo- 
geneous atmospheres utilizing a Henyey-Greenstein phase 
function or a Modified Henyey-Greenstein phase function [16], 
and a given portion of Rayleigh scattering. It was necessary 
to modify the routine so that it adequately represented the 
geometry of a finite cloud nd so that any arbitrary phase 
function could be adapted to it. These changes are summarized 
in the following sub-sections. 

1. The Model 

Figure 4 represents the model used in this simulation. 

The model is identical to that of Ref. 6 with respect to photon 
path ration and accountability. The figure depicts a 
cloud on whose left boundary the photons are incident at the 
center of the accountability shells, and at whose right 
boundary time and spatial information is tabulated. Figure 5 
expands the cloud and lists the names of the parameters used 
in the simulation.  KSCAl and KSCA2 represent the scattering 


Mefficients in inverse kilometers of the inside and outside 
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media respectively. Likewise, KABSI and KABS2 represent the 
absorption coefficients of the two media. RPTI and RPT2 
represent the ratio of particulate to total scatter and Gl 

and G2 represent the Henyey-Greenstein phase function para- 
meters when used in this simulation. THICK is the physical 
thickness of the cloud in kilometers. Appendix C gives details 
on the actual modification of the routine. 

At each crossing of a photon from inside of the cloud 
to beyond the cloud, the total distance traveled is determined 
and tabulated in time of arrival bins. Time spread and delay 
information is presented using these bins. A running summa- 
tion and counter are used to calculate the mean and standard 
deviation of time required to reach the cloud boundary. 

One other change from the original routine is conversion 
of dimensions from units of extinction lengths to units of 
kilometers. Only minor changes in program logic were required 
to make this change. 

Typical values of atmospheric parameters were used in 
most simulations. A summary of the parameters used is given 
in Table I. All diagrams of spatial and temporal character 
include parameters used in the specific case. 

The Monte Carlo results of Bucher [1] were very useful 
as guidelines for determining correct operation of the model 
created. Comparisons of the two models appear later in this 


section. 
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TABLE I 


Scattering Coefficients used in the Simulation 


1 


KSCA(km `) Albedo 
Clear Atmosphere „291 1.0 
Light Haze . 90 120 
NOSC Fog 6: 57 Lc 
Water Cloud C.2 11.18 1.0 


2. Phase Functions 
As mentioned in the statement of the problem, three 
phase functions of different peakedness yet equal in (cost) 
were needed to investigate the problem. The following phase 
functions were employed: 
a. Henyey-Greenstein Phase Function 
Figure 1 depicts the functional dependence of 
scatter on angle off direction of incidence, 0. The closed 
form expression for the Henyey-Greenstein function has been 
stated earlier in Equation 8. A G value of .83 was selected 
so as to match other phase functions used in the simulation. 
Ihe peak of the phase function is only five units. Reference 
6 explains its use in the Monte Carlo routine. 
b. Numerical Data Input for Arbitrary Phase Functions 
The other two phase functions required development 
Of a method for inputting phase functions consisting only of 
data pairs and no closed functional form. They are that 


representing Water Cloud C.2 [21] at .53 microns and that of 
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NOSC Fog. The first was generated by use of the program of 
Appendix B, the second by a similar calculation using an 
experimentally obtained water particle distribution. Appendix 
A describes the details of generating a random scatter angle 
weighted by these arbitrary phase functions defined by data 
pairs. 

c. Rayleigh Phase Function 

This routine can also simulate a fractional amount 
of the total scatter as Rayleigh scatter. In this work the use 
of this capacity was minimal because in most cases water clouds 
of predominately large particles were simulated. 
3. Automation of Contour Plotting 

Contours of equal relative photon flux, as described 
in Ref. 6, were used in this work for spatial characterization. 
A method of machine computation was adapted to the DRLITE 
routine. This allowed the computer to perform the tedious 
interpolation made before by pocket calculator. Although 
useful in some instances, the procedure proved quite computer 
time consumptive. Its usefulness was very pronounced in comparing 
the Monte Carlo routine with Gordon [9] as will be discussed 
later. Appendix D gives details of the adaption of this method 


ES the routine. 


B. ANALYTICAL METHODS 
As explained in the previous section on theoretical pre- 
liminaries, analytical methods ordinarily lead to solutions 


Which are applicable only under certain conditions. This is so 
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because in order to reach a usable solution, approximations 
have to be made. Nevertheless, in the regions of their 
applicability, analytical solutions are excellent sources of 
the information needed to verify results generated by other 
more universally applicable methods. In this work two analyti- 
cal treatments were employed for just that reason.  Gordon's 
[9] paper on practical approaches to underwater scattering 
Was used to verify the spatial output of the Monte Carlo routine 
created and Stott's [26] closed form expression for time 
Spread was used to verify the temporal output of the program. 
The two treatments are summarized in the following paragraphs. 
1. Effective Attenuation Coefficient (EAC) Method 

Gordon's paper [9] presents the concept of an effective 
attenuation coefficient which considerably simplifies under- 
Water multiple scattering calculations. Concise expressions 
for the total flux through an aperture and the beam-spread 
function are defined in terms of the EAC. The derivation of 
these expressions is included for easy reference as Appendix E 
as is a documentation and description of the program adapting 
them to machine computation. It is interesting to note that 
Gordon's treatment is actually a truncated version of the 
exact analytical solution of the multiple scattering problem 
using the small angle approximation provided by him [28]. The 
truncation, in effect, considers the scattering phase function 
as a forward scattering delta function, thus does not consider 
higher order terms of the exact solution. It is not surprising 


then that the result is only applicable under certain conditions. 
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Fortunately these conditions are found often in the underwater 
optical world. The comparison of this method to the Monte 
Carlo method is found in the next section. It should be kept 
in mind that Gordon's treatment, even in the truncated manner, 
closely agrees with experimental data [29]. The reader is 
referred to Appendix E for further details on the EAC method. 
2. Closed form Tıme Spread Expression 

Stott's treatment [26] of time spread in a multiple 
scattering region is short, concise and seems to agree well 
with existing time spread data. It was therefore chosen as 
an additional means for verification of time spread results 
given by the Monte Carlo method. The details in deriving 
Stott's closed form solution for time spread are located in 
Appendix 6 for easy reference. The result is the average 
additional multipath time required for a photon to traverse 


a scattering medium with physical thickness z, optical thick- 


Bess t, albedo My and RMS scatter angle y,. The time spread 
is 
L {с [G*z as ty 501 5-1]-1 (12) 
C VootYo 


This result is compared to Monte Carlo results in later sec- 


mons. 


©. COMPARISON OF METHODS EMPLOYED 
A problem usually encountered in mathematically modeling 


any physical situation is to ensure at all times that the model 
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employed is turning out reasonable answers, To verify that 
this was in fact the case, frequent cross-checks were made 
between the models. The proof of the theory follows, of 
course, in agreement between the models and/or agreement with 
experimental observations. This section compares the results 
of the Monte Carlo model with the analytical models as well 
as with the results of Bucher's [1] independently created 
Monte Carlo model. Once confidence is established, conclu- 
sions can be drawn about the regions of effectiveness of each 
model. Other more routine checks for error are noted in 
Appendix F. 
l. Monte Carlo vs. Analytical 
a. Spatial Characterization 

The EAC Method derived and verified in Appendix E 
is compared to the Monte Carlo method in this section. First, 
the method of spatial characterization adapted from Ref. 6 is 
briefly explained. 

(1) Description of Relative Flux Contours. At 
each angular bin of each reference shell the relative flux is 


given by 


# Photons Passing Through Bin 


Total Number Ejected x Area of Bin (13) 


Relative Flux = 


The negative logarithm of this quantity is calculated for each 
bin. This array of values is interpolated as described in 
Appendix D,and contours of equal relative photon flux are 


drawn. Likewise, the negative logarithm of the beam spread 
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function evaluated by the EAC method is evaluated at a similar 
array of spatial positions and equal flux contours are drawn. 
The two sets of contour lines are compared in the following 
paragraphs. 

(2) Comparison. Flux contours calculated using 
Henyey-Greenstein G parameters of .70 and .95 for each albedo 
of .80 and .90 were drawn at integer values of the negative 
logarithm. Figures 6 through 9 show these comparisons, The 
dotted flux contours are those of the EAC method from zero to 
45 degrees and the solid contours are those of the Monte Carlo 
method. The term S on the graph is the scattering coefficient 
as well as the albedo since the extinction coefficient is 
unity in each case. Dimensions are in kilometers but.in this 
case one kilometer 2 equivalent to one extinction length. 

In general, agreement is good between the 
methods for highly peaked, G - .95, phase functions at dis- 
tances less than 10 extinction lengths and angles varying from 
one to 20 degrees. This comparison agrees with results given 
in Gordon [9]. It should be noted that angular agreement is 
nearly independent of albedo whereas range agreement is better 
when large absorption is present. 

(3) Conclusions. Inside of a few degrees at moder- 
ate to large extinction lengths the EAC method becomes 
undependable in most cases. The Monte Carlo and EAC methods 
agree well in regions where agreement should be expected. The 


Monte Carlo method has the advantage over the EAC method when 
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considering atmospheric scattering because in this case 
absorption is small and large optical thicknesses are encoun- 
tered often. Confidence can be atea in the Monte Carlo 
routine for spatial characterization of light in a multiple 
scattering region. 

b. Temporal Characterization 

Stotts's [26] closed torm expression for time spread 
which is summarized in Appendix G is compared to time spread 
found using the Monte Carlo method in this section. The 
technique used in tabulating time spread information in the 
Monte Carlo model was explained earlier in the section titled 
The Model", 

(1) Comparison. Multipath time spread is defined 
as the difference between the average transit time incurred 
from multiple scattering and the normal transit time in the 
absence of scattering. The expression for time spread derived 
by Stotts is given in Equation 12. This equation is plotted 
in Figures 10 through 14 with an albedo of one, an RMS scatter 
angle of .65 and optical thickness, t, ranging from zero to 
NO for physical thickness values, z, of .5, 1.0, 2.0, 3.0 and 
4.0 kilometers. Monte Carlo data points are marked with X's 
using various t values with other parameters the same. Also 
plotted are Bucher's [1] data which is explained later. In 
general, agreement is within one microsecond for optical thick- 
nesses less than 10 and within five microseconds for optical 


thicknesses between 10 and 60. 
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(2) Conclusions. The Monte Carlo and closed form 


methods agree well in regions where agreement should be expec- 


|ted. Because of this agreement, confidence can be placed in 


the Monte Carlo routine for temporal characterization. The 
Monte Carlo routine has the advantage over the closed form 
expression in that different phase functions can be simulated, 
thus allowing study of the effects of the details of the phase 


functions on time spread of light in a scattering medium. 


Because time spread is expressed quite well by the closed form 





expression, little dependence on these details is expected. 
The dependence is discussed in later sections. 
2. This Monte Carlo Model vs. Bucher's Monte Carlo Model [1] 
a. Spatial Characterization 
The Monte Carlo method is tested for agreement 
with Bucher's [1] Monte Carlo results. Because both methods 
are built with nearly the same geometry, nearly identical 
results are expected. First the graphic representation of 


data used by Bucher needs explanation. 


(1) Description of Scaling Method. Bucher [1] 


defines for a scattering medium the diffusion distance 


D 


Pa © Zosis ° 


(14) 


Where D is the mean free path between collisions and (cos9» is 


the average cosine of the scattering angle as before. Using 
this formula, the effective scattering thickness ta is given by 


the relation 
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таб” t(1-(cos0») (15) 


where T is the physical thickness of the cloud and T is the 
| usual optical thickness. The spatial results are then pre- 
sented by graphically displaying the average displacement 


(r) in relation to the diffusion thickness, D as a function 


ar 


of effective scattering thickness Ta: Also presented is the 


critical transverse displacement, Po? in relation to Da inside 


— 


which 50 per cent of the photons exit the cloud at effective 
scattering thickness, Tae The reader is referred to Bucher 
[1] for further details. 
(2) Comparison. Monte Carlo calculations were 

| made and interpreted in terms of the units described above. 
| Figures 15 and 16 show the results at various effective 
Scattering thicknesses using an albedo of unity and a Henyey- 
Greenstein function for which G = .83. The X's represent data 
collected here and the continuous line is a best fit curve to 
Bucher's data using the same parameters. As can be seen, 
| agreement is excellent in both cases. 
(3) Conclusion. Spatial characterization by this Monte 
Carlo routine is consistent with a previously developed routine. 
| This creates confidence in the routine. It can now be used to 
| Study the effect phase function details have on the spatial char- 
Heter of light. 
b. Temporal Characterization 

The time spread data of the Monte Carlo routine is 
compared to Bucher [1] in this section. Two methods of 


— — 
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Figure 15 
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Figure 16 


This Monte Carlo vs. Bucher's 
Monte Carlo [1] Spatial Spread Comparison 
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presentation are used. One uses time bins, the other uses 


Bucher's best fit equation for time spread. 
(1) Comparison. Bucher's equation for multipath 
time spread as defined in an earlier section is, 


95 1 
ys 627, ): 


_ 10 E 
L(t е = 3.91 x 10 T(t 3) (I6) 


d 


| where L(t.) is the time spread in seconds, T is the physical 
thickness of the cloud in meters, Ta is the effective scatter- 
ing thickness of the cloud and c is the speed of lisht in 
meters per second. This function is plotted alongside Monte 
Carlo data and Stott's closed form equation in Figures 10 
throush 14 for an albedo of unity, a Henyey-Greenstein phase 
function with G = .83, and a physical thickness of .5, 1.0, 
2.0, 3.0 and 4.0 kilometers. 

Agreement is well within one microsecond for 
nearly the entire range which is well within statistical 
error. In these trials only 50 - 500 photon histories were 
tabulated which seems to be an indication that collection of 
time spread information does not require unwieldy amounts of 
computer time. 

Another comparison of time spread information 
is graphically displayed in Figures 17 and 18. The number of 

photons arriving in each time bin is normalized to the maximum 
| number in any bin and plotted as a function of time in micro- 
| Seconds. Bucher's data is superimposed by a dotted line on 

the time bin data presented. Figure 17 plots results for an 
Optical thickness of 15 and Figure 18 for an optical thickness 
of 30. Other parameters are listed on the graphs. Agreement 
Once again is very good between the methods. 
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(2) Conclusion. Once again, confidence is ensured 
in temporal characterization of time spread and pulse spread 
by the Monte Carlo routine. The routine can now be used in 


investigating time spread dependence on phase funtion. 
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V. RESULTS 


A. SPATIAL SPREAD DEPENDENCE ON PHASE FUNCTION 

In this section five graphs are used to display the 
dependence of the spatial spread of light on phase function 
when traversing a scattering medium. All five diagrams 


show the relative flux contours described earlier in this 


thesis and in Ref. 6. Each figure lists the parameters and 
the phase function used in the simulation. Each axis is in 
units of kilometers. Incidence of photons is at the origin 


and to the right in all cases. The left edge of the cloud, 
if it exists, is along the axis perpendicular to direction 
of incidence and the right edge of the cloud is indicated by 
a dotted line if it is within the boundaries of the graph. 
For easy reference, Table II tabulates the parameters and 
phase function used in Figures 19 through 36. 

The phase function of Water Cloud C.2 is used in Figure 19 
to show the peakedness of the relative flux contours out to 
at least 20 extinction lengths when the albedo is only .90. 

A "toe" appears on the contours at small angles due to the 
, high degree of forward scattering of Water Cloud C.2. 

Figures 20 and 21 compare the flux contours due to NOSC 
Fog and Henyey-Greenstein using G = .83, respectively. The 
flux contour for NOSC Fog is more forward peaked than that of 
the Henyey-Greenstein function for distances less than three 


Kilometers or 19 extinction lengths but flux contours are very 
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Similar at distances greater than three kilometers. Outside 
of 20 extinction lengths the contours and therefore the 
spatial character of the light they represent become indepen- 
dent of the phase function used. The contours of these two 
graphs are much less forward peaked than that of Figure 19 
because the albedo in the former is unity. 

Figures 22 and 23 show similar results when comparing 
results when Water Cloud C.2 [21] and Henyey-Greenstein were 
used, respectively. In all figures the cloud was extended to 
infinity at the right of the origin, explaining the discontinu- 


ity in contour slope. 


B. TEMPORAL SPREAD DEPENDENCE ON PHASE FUNCTION 

In this section nine histograms and one table are used to 
display the dependence of time spread of light on phase func- 
tion when traversing a scattering medium. All nine diagrams 
list the parameters in the same fashion as in the last section. 
The essential parameters can also be found in tabulated form 
in Table II. The information is plotted as the number of 
photons arriving in any bin normalized to the bin of most 
frequent arrival versus the time spread incurred. Note that 
the bins are logarithmic and therefore represent much less 
Span in time to the left than to the right. 

Figures 24, 25 and 26 compare time spread for NOSC Fog, 
Water Cloud C.2 and Henyey-Greenstein, respectively in a cloud 
ene kilometer thick of optical thickness four. It can be seen 


by comparing the graphs that at this optical thickness the 
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higher peaked phase functions have more photons arriving at 
earlier times. 

Figures 27, 28 and 29 similarly compare time data for a 
cloud with a physical thickness of one kilometer and an optical 
thickness of eight. The time spread diagrams begin to look 
very much alike but still vary slightly with peakedness. At 
this thickness the average transit time is nearly identical 
for the three cases but there are still photons of the higher 


peaked phase functions which traverse the cloud with very 


little time delay. This indicates some dependence on phase 


function for this thickness. 

Figures 30, 31 and 32 conclude the graphical time spread 
presentation with a cloud one kilometer thick with an optical 
thickness of 15. There is no noticeable difference in the 
graphs that can not be attributed to the statistical deviation 
expected in a Monte Carlo routine. Table III summarizes the 
time spread data. 

From these results it can be said that the time spread 
becomes independent of the details of the phase function at 


optical thicknesses greater than 15. 


C. REGION OF TRANSITION FROM FORWARD TO MULTIPLE SCATTER 

The region of transition from forward to multiple scatter 
can be defined as that region for which (1) distances of lesser 
magnitude show time and spatial spread depending markedly on the 
phase function of the scattering particles and (2) distances of 


greater magnitude show that time and spatial spread do not depend 


— M 


on the phase function of the scattering particles. 
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TABLE Phi 


Summary of Time Spreading 
by One Kilometer Thick Cloud 
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From the previous two sections the region of transition 
from forward to multiple scatter is between about 10 - 20 


optical thicknesses. This agrees favorably with observations 


' made by other authors [40]. 


D. EFFECT ON SPATIAL CHARACTER OF LIGHT WHEN PASSING THROUGH 
A CLOUD 


In this section four graphs are used to display the effect 
of finite thickness clouds on the spatial spread of light. 

All four diagrams use the relative flux contours and the list 
of parameters defined in earlier sections of this report. 

Figure 33 shows the relative flux contours of an infinite 
very clear atmosphere. Figure 34 shows the effect when a cloud 
of optical thickness about nine is placed in the light's path. 
Deformation of the contours is very obvious. 

Likewise, Figures 35 and 36 show the effect of a .5 kilo- 
meter Water Cloud C.2 when placed in the path of light in a 
clear atmosphere. The sharp peaks on the flux contours are 
rounded and as expected, light is significantly redirected 
when passing through the cloud. Clouds of large optical thick- 
ness would eventually cause the Lambertian irradiance well 


known in the theory of light propagation. 
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VI. DISCUSSION 


In this thesis Monte Carlo methods and simple analytical 
methods were used to characterize light transmission through 
a multiple scattering medium. Complicated numerical methods 
were not used nor were highly mathematical models. This type 
of approach remains for future investigations. It is antici- 
pated that results drawn from present methods will be confirmed. 

Much information on the multiple scattering problem re- 
mains to be gathered. The present Monte Carlo routine can be 
adapted to numerous useful geometries and atmospheres with 
gradients and time variation of physical parameters.  Propaga- 
tion over and around solid bodies such as mountains and the 
earth itself can be modeled.  Layered media such as cloud, 


air and water interfaces can also be simulated by the routine. 
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Vit. CONCLUSIONS 


A computer routine based on a Monte Carlo model was de- 
veloped to simulate light propagation through a scattering 
absorbing medium. The routine can simulate a plane parallel 
cloud of scatterers with the desired parameters located within 
a medium with another set of parameters. The phase functions 
of each medium can be arbitrarily defined using a set of data 
pairs or can be approximated by closed form expressions. The 
data generated by the routine has been checked for accuracy 
against other theories using analytical methods. Each compari- 
son shows adequate agreement between theories where agreement 
can be expected. 

The routine will automatically plot spatial information 
necessary in characterizing light transfer through the medium 
by use of contours of equal photon flux. It will also generate 
histograms depicting time spread information for light passing 
throush finite clouds. 

Using the Monte Carlo routine created and inputting phase 
functions of different peakedness, it has been found that both 
spatial and temporal spread are independent of the details of 
the phase function for thicknesses greater than 15 extinction 
lengths. The region of transition from forward scatter to 
multiple scatter is between 10 - 20 extinction lengths. 

The routine has also been used to study the effect finite 


thickness clouds have on the spatial character of light. 
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APPENDIX A 

I. GENERATION OF A RANDOM SCATTER/ANGLE 

WEIGHTED BY AN ARBITRARY PHASE FUNCTION 
me, STATEMENT OF THE PROBLEM 

Reference 6 outlines the theory and calculations necessary 

to generate random numbers weighted in accordance with functions 
that represent characteristics of a scattering medium. At each 
collision new scatter angles and a distance were generated by 
closed form expressions which weight a uniformly distributed 
random number. However, in the more general case of a poly- 
dispersion, the computed phase function could not be 
represented adequately in closed form. This made it impossible 
to invert the function enabling analytical generation of a 
weighted scatter distribution. It was necessary to create a 
method for generation of a random theta weighted by an arbitrary 


phase function. 


B. METHOD USED IN WEIGHTING THETA 

Reference 21 and Appendix B explain in great detail a 
method for calculating the representative phase function of a 
polydispersion given the particle size distribution and com- 
position. Using this computer adaptation of Mie theory, values 
of the averaged normalized phase function were generated at 
selected angles of scatter. Given this data and a random 
number, R, weighted uniformly over the closed interval [0, 1], 


the problem is to solve the following equation for E 
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(1) 


where P(8) is the averaged normalized phase function and 0 is 


the variable of integration. Notice that, 


CD 
и 


О when R 


O, 


Ө = т when R = 1, 


as you would expect on the closed interval.  P(8) is not a 
continuously defined function over the interval [0, m] in this 
case so the integrals are represented numerically by the 
trapezoidal rule. Figure 37 diagrams the basic procedure used. 
N values of P(8) are selected so as to closely approximate the 
phase function. Of course, more values are selected near the 
small scatter angles to adequately describe the sharp peak. 
The N values are multiplied by the sine of the respective 
Scatter angle (P(0) includes sin 0 implicitly in Figure 37). 
N-1 trapezoids are established using these N values. The area 
of each trapezoid is divided by the total area of all trape- 


Ezoids thus establishing a weight for the respective interval. 


"- (P:sin0;*P.:,,Ssin0;,,)(0;,,-0,) em. 


i N-1 
Y (P;sin0i*P:,,Sin0:,12(9;,41-01) 
I l 


This of course requires that 


W s (4) 
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GENERATION OF A RANDOM THETA 
WEIGHTED FOR AN ARBITRARY 
PHASE FUNCTION 


Points generated by 
Computerized Mie Theory 


N-f trapezoids are constructed 


P (0) from N data pairs 
Ва умт 
| ie wo 
ip 
i a Жыр — 
Pott (Or Pan 
Area of 
- W 
W, Me ч Trapezoid N-1 
number i 
9, Ө 


EA 4 
R= Area between lines 

R = Area of portion ү 

between these lines 

under curve 








Figure 37 


Diagram of Random Generation of Theta 
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With the weight of each panel known, the uniform random 


number R is used to solve for the first value of n such that, 


n 
2, W. > R (5) 
i=1 


is satisfied. The random Op required is somewhere between 9 


and a Within the panel, 9 is weighted linearly in a fashion 


+] ` 
as explained in Ref. 6. For a given panel, E 9:41? Р. апа 
Р. т are known. From these the slope and intercept are estab- 
lsihed. 
P.,,-P. 
В = — A. = P.-B.6. (6) 
1 sse 1 Jes L 
er] 


The panel is normalized as follows: 


9i+1 
NORM / (A,*B,0)d8 = 1 (7) 
тү 
Ө. 
1 
so that, 
1 
Anton ODE See] -9; ) 
А11 that remains is to solve the following for 9 
n 
Ra a 2 мокм/ ХА, +В,9)909 (9) 
(y. LR. y.) m 
E 1 
1=1 1=1 
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After a few algebraic steps, the solution is 


1 





See Been (10) 
R B. De Be - B.«0 
1. i 1 ai 
where, 
B.6.,° B,0,? 
C: = 2, (8,0452) | 2 +А.0. CIT) 
This On has the desired properties. The FORTRAN coding of this 


method is found in the Subroutine RANTH. 
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APPENDIX B 

I. DESCRIPTION ANDODOCUMENTATIONGOE PROGRAM 

TO ADAPT, MIE THEORY TO MACHINE COMPUTATION 
A. INTRODUCTION 

Mie theory has been adapted to machine computation on many 

occasions.  Deirmendjian [21] provides an excellent guideline 
for using Mie theory on spherical polydispersions. Using 
Deirmendjian's guideline, a computer routine was created to 
generate the volume scattering and extinction cross sections 
and the corresponding elements of the normalized scattering 
Matrix for a polydispersion where the number of particles per 


unit volume, per unit radius is given by, 


OUTILS 
n(r) - ar e ES o<r<w СШ) 
where r is the particle radius. The four constants a, a, b 
and y are positive and real and «a is an integer. They are not 


independent of each other, and are related to quantities in the 
particle frequency distribution. The radius which is most 
Frequently found in the particle distribution i5 ГЕ апа N is the 
total number of particles per unit volume. Both N and r can 

be found by experimental measurement. In terms of N and D. the 


constants of the distribution are found using: 


S o -br' a atl 
Nasa re dr = — əƏ> T(—) 2) 
О 
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Y | (3) 


and by choice of a and y to best fit the experimental distri- 


Bution function. T is the usual gamma function. 


B. EQUATIONS TO BE CALCULATED 

The equations used to generate scattering data fora 
polydispersive system are extensions of those used in a mono- 
dispersion. The distribution has been created in a manner 


requiring that 


2 
N Ji n(r)dr (4) 


р 
JE 
where n(r) is a continuous and integrable function within the 


range and represents the partial concentration per unit volume 
per unit increment of radius r. The values of interest are 
volume scattering cross sections and the corresponding elements 
of the normalized scattering matrix [18, 21]. These values 


can be computed using 


sca 


со 
В [PX п хх) =š nef х?п(х) Оѕса(х) ах (5) 
O 


and 


00 
Вех А» n(x) J nef X а(х) ОБ оох (6) 


Р; (0) = rises] nex) і. (0) ах J 1.2 (7) 
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where 


Q (x) = 


Sca 


Q (x) = 


ext 


i.(x) = 
J 


scattering efficiency factor 
extinction efficiency factor 


dimensionless intensity parameters 


Each of the above terms will be examined in detail in following 


sections. 


C. CALCULATION OF SCATTERING EFFICIENCY FACTORS AND DIMENSION- 


LESS INTENSITY PARAMETERS 

To perform numerical integration, the 
evaluated at many different values within 
In this case the scattering and 


gration. 


factor must be evaluated for numerous Mie 


integrand must be 
the range of inte- 
extinction efficiency 


Size parameters x. 


Scattering efficiency is a name given to the ratio of total 


Scattering cross section to geometrical cross section. A 


Similar definition is used to define extinction efficiency. 


For any particle size and composition, Mie theory gives the 


total scattering cross section as [18, 21], 


ô cca MX? : J (A, AS FA, Af) do (8) 
Q 
where 
kA, = S, (m,x,0) - X D Gub? (9) 
2 2n+1 
KA, = S,„(m,x,0) = 2 ACHE соата “пл” (TD) 
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using the complex index of refractıon of the particle m, the 


Mie coefficients aj,» b, and angular coefficients Ta? T each 


II ° 


of which is described in detail in later sub-sections. Using 


these equations, the scattering efficiency factors are [18, 21], 


CD 
2 
Qa oí 0053 7 GF 2 (2n+1) (la, 1*+]b,12)2, (11) 
x n=l 
2 00 
Qayt (MX) = = 24 (2n+1)Refa +b }. (12) 


The dimensionless intensity parameters are given by the 


expressions, 


e = 2 RH = de 
i, (x,m,6) AK AAF S191 

(13) 
: = k? # = zs 


A S, and S, are the same as above. The procedure 


1 2 2 EE 2 


used in calculating E b 


where A 
nomm and T is described next. 
l. Computation of та апа т 

The two angular functions e and T, are required before 
the dimensionless intensity parameters can be calculated. t. 
and m are functions of u = cos8 only and as the subscript 
implies, there are actually many different angular functions. 
Both functions are defined in terms of Legendre polynomials [18] 


and their derivatives. The following recursion relations are 


used in calculating u and T n» 
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_ 2n-1 1 
m 08) = cos T mq 1 9? 7 7 aa (14) 


r (8) = [cose perp] - [(2n-1)sin*8 m _,(6)] 
+ T2290? (15) 
with 
т (60) = 0 т (0) = 0 
о о 
7, (8) = 1 Tt (8) = сос Ө (16) 
O 
7 (8) = 3 cos® Tt, (8) = 3 cos 26 
2. Computation of a_ and b 
n n 


The Mie coefficients a, and b. have a complex depen- 
dence on the index of refraction of the particle, and the 
surrounding medium and also depend on the size parameter, x. 
There exists many different expressions for the Mie coeffici- 
ents in terms of known mathematical functions. The form used 
in this program is described in terms of spherical Bessel 
functions each of which is described in later sections. The 


Belatıons for an and b, are 


ү H уг a 00723, Cy] x3, (x)-ey] ¿[xi 100-7 -nj (x) | en 
D [yá 1 9-23, 69] O | x) | 
BE SEN 1 (y)-nj КЛЕИ (x)-uyj sty) [xd 1“Х)- -nj К] 
j Alm 1552-7 nj пб хва (x)-uyj a [nz Fan G0 | (18) 
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where is the absolute value of the complex index of refraction 
outside the sphere, € is the absolute value of the complex 
index of refraction inside the sphere, x = 2mr/iX and y - mx. 


The functions E VOS l J,71Y, are the spherical Bessel and 


n 


Neuman functions. A procedure for calculating their values is 


outlined in the following sections. 


2 


3. Computation of > Ah, 


ү? 


The argument of these Bessel related functions is 


complex in some cases, so much care is required in their calcu- 


lation. Each function is calculated using recursion relations 
and the lowest order functions as follows. For jn? 
(2n*1)j, (z) 
Ја+102) 2 77 5 13-162) B 
: _ Sinz 
Jg C2) аба 
: . Sinz _ cosz 
ji (2) a (20) 
3562) = ( е Z)sinz E cosz 
= 2 


likewise for s 


Е (2п+1)уһ(2) 


У. 2 = y, 412 (21) 
y "= 

yy (2) = “2052 . Binz (22) 
Y, (Z) = pe z) cosz - —Ñ sinz 


EDI 





and these two yeild hos 
2 BIER: _- 
h (2) = jj Cz? іу (2) (23) 
Because the argument of the sine and cosine in evaluating the 


spherical Bessel function is complex, the following relations 


are needed, 


z sog d T5 
cosz = cos(a)cosh(b) - isin(a)sinh(b) (24) 
Sinz = sin(a)cosh(b) + icos(a)sinh(b) 


where cosh and sinh are the usual hyperbolic trigonometric 


Functions. 


D. DESCRIPTION AND DOCUMENTATION OF PROGRAM 

The program adapting Mie theory to machine computation is 
composed of the MAIN routine and numerous subroutines described 
in the following sub-sections. 

1. MAIN Routine 

The MAIN routine handles the input and output funetzons 

necessary in program operation. The input parameters include 
indices of refraction inside and outside of the spheres, con- 
stants of the particle distribution and the most frequently 
occurring value of Mie size parameter in the distribution, x. 
Through Xo» the wavelength of the incident light is entered 


because p is known for any desired particle distribution. 
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Other necessary inputs are the number of scatter angles desired 
in the output scattering matrix elements ZU the smallest 
value of x, the increment in x to be used for numerical inte- 
gration and the odd number of x values at which the integrand 

is to be evaluated. The MAIN program accepts the input values 
and uses the distribution parameters to assign each designated 
X value a weight. Because the particle distribution is normali- 
zed, the sum of all weights assigned is unity. 

Simpson's 1/3 rule is used to evaluate integrals 5, 6 
and 7. At each value of x the scattering and extinction 
efficiency factors of equations 11 and 12 are calculated through 
use of a subroutine MIEM. Similarly, at each desired scatter- 
ing angle the dimensionless intensity parameters of 13 are 
calculated through MIEM. Printout can be called at each value 
of x if so desired. The MAIN program cumulatively sums the areas 
of an even number of panels to get the desired results which are 
then printed. 

2. MIEM Subroutine 

MIEM is a subroutine composed to calculate the scatter- 
ing and extinction efficiency factors in the integrand of 
equations 5 and 6 and the dimensionless intensity parameters of 
equation 13. Required inputs of MIEM, transferred to it by 
MAIN, are the index of refraction of the particle (henceforth 
the index of refraction of the outside medium will be unity), 
the size parameter x, and the number of scattering angles at 


which calculation of the intensity parameters is desired. MIEM 
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uses equations ll and 12 to calculate the efficiency factors 
and equation 13 to calculate the intensity parameters. 
Because equations 9, 10, 11 and 12 require summation of series, 
each element of the series must be evaluated in turn. The 
terms involved in finding the efficiency factors are an and 
b. which are calculated for each n * 15 2, 3, 5... úntil 
the next term of the series is adequately small or the total 
number of terms exceeds 120. After each calculation of a, and 
b. its contribution to "MSS and Q est X? is added to the pre- 
vious total.  MIEM uses function subroutines JN, FN and HN to 
calculate each of the a, 5 and b 's. Arrays JX, FX and HX are 
used to store values of the corresponding spherical Bessel 
functions during each recursion step. 

After all the parameters a, and b. are calculated, 
MIEM turns them over to ANGLE to complete the remaining dimension- 
less intensity parameter calculations. 

3. ANGLE Subroutine 

ANGLE is a subroutine called by MIEM to calculate the 
dimensionless intensity parameters at each desired angle of 
scatter. ANGLE requires as inputs the Mie scattering coeffici- 
ents generated in MIEM and the number of scatter angles desired. 
When called, ANGLE uses equations 9 and 10 to find S, and 5, at 
each scatter angle. To do so ANGLE evaluates en and Ta using 
recursion relations 14 and 15 with initial order functions, 16. 
At each n, the corresponding Mie coefficients an and b. are used 


with E and T and their contributions to 5+ апа 5% are added. 
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The process is repeated for each scatter angle and the result- 
ing arrays are used to evaluate equation 13 at each scatter 
angle. The resulting dimensionless intensity parameters are 
then returned to MIEM and in turn to MAIN for integration using 
equation 7. 
4. JN, FN, HN, CCOS, CSIN Complex Functions 

These functions are called by MIEM in evaluation of 
the Mie coefficients a, and bo: JN, FN and HN contain logic 
to perform the recursion operation of equations 19 through 23. 
CCOS and CSIN are complex trigonometric functions drawn upon 


as needed by JN, FN and HN. 


m CPU TIME CONSIDERATIONS 

CPU time depends largely on the wavelength to particle size 
ratio. This, of course, varies with the particle size distribu- 
tion used. Thus, if the distribution includes many particles 
of large size compared to the wavelength the CPU time is great 
and vice versa. Typical CPU time requirements were on the order 
of 15 to 30 minutes for wavelengths of .53 to .28 microns, with 
about 300 x values in a distribution of Water Cloud C.2 of 
Deirmendjian [21]. The dependence on particle size to wave- 
length ratio is due to the fact that many terms are required for 
convergence of series for large values of x. As expected, CPU 
time goes up guite linearly with the number of x values used as 
integration points. There is little dependence on the number of 


scatter angles required. 
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APPENDIX C 
L. A SIMULATION TOF AxCLOUDIN 
MONTE CARLO ROUTINE LITE 

A. INTRODUCTION 

Reference 6 describes Monte Carlo simulation of light propa- 
gation through an infinite, homogeneous atmosphere. Many 
problems can be sufficently investigated using this model, but 
one of the advantages of a Monte Carlo model is its relative 
easein adaptation toirregular geometries and inhomogeneous 
atmospheres. This Appendix gives one method for simulation of 
a cloud with plane parallel homogeneous medium. All macro- 
scopic properties are the same everywhere inside of the cloud 
and another set of macroscopic parameters are the same every- 
where outside of the cloud. This model circumvents the entire 
problem of describing and locating boundaries and inhomogeneities 
in real clouds. Figure 5 depicts the general structure of the 
cloud model giving the names of various parameters of the com- 
puter simulation. A point source of light is incident normal 
to the left edge of the cloud and its path is randomly generated 
until it exits the outermost sphere of the reference volume. 
Reference 6 describes the random path generation and accounta- 
bility also used in this model so the terminology of that 


reference is used here for continuity whenever possible. 
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B. MODELING CLOUD BOUNDARIES AND PROPERTIES 

The sets of parameters needed in defining two different 
Scattering media are implemented by use of storage arrays. A 
binary logic code is used that switches whenever a photon 
crosses a boundary. Each array of characteristic parameters 
has two columns, one for inside the cloud and one for outside 
the cloud. The logic switch determines at each collision 
from which column the parameter is to be drawn. 

Upon crossing a boundary, the photon is stopped and a new 
distance is randomly generated in accordance with the para- 
meters of the newly entered medium. The method used for 
determination of whether or not a boundary was crossed is 


described in the following section. 


C. DETERMINATION OF BOUNDARY CROSSING 
The following equation is used to determine at each colli- 
Sion the angle between the original direction of incidence and 


the present position vector, R, 


S B R-(R-R, ) 
= cos — 

r 
where r is the distance between the collision point and the 
point of incidence on the cloud. Ry 1S an orientation vector 
from the point of collision to the (0, 0, 1) point of the fixed 
coordinate system. The r and @ values are computed at each 


collision which allows, due to cylindrical symmetry, calculation 


of the projected distance along the direction of incidence. 


10.2 





This distance is 


Projected distance = r cos 0 = DIST, 


which is compared to the thickness of the cloud. The medium 


in which the photon exists is determined by, 


DIST ow < 0 Behind 


UE DIST « THICKNESS(Inside 
new 


DIST" > THICKNESS| Beyond 


ew 


The location of the collision relative to the cloud is compared 
to the location of the previous collision and a boundary 
crossing is found if it has occurred between collisions. Logic 
was created to then stop the photon at the boundary and project 
it along the same path using new scattering parameters. There 
are nine different combinations of old and new positions, three 
possibilities for the old position and three possibilities for 
the new position. Each specific situation is investigated at 
each collision and upon boundary crossing the correct stopping 


formula is applied. 
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APPENDIX D 
I. AUTOMATION OF RELATIVE FLUX CONTOUR 
PLOTTING IN DRLITE ROUTSNE 

A. INTRODUCTION 

One very important output of the Monte Carlo routine is the 
relative flux per unit area at numerous grid points within the 
reference spheres with respect to the original energy projected 
by the point source. Knowing the relative flux at numerous 
grid points allows generation of equal photon flux contours. 
Calculation of the contours was a tedious job requiring many 
hours on a programmable pocket calculator and at the drawing 
board. Accuracy of the contours suffered (not to mention the 
one doing the work) therefore a method was created to produce 


many of the contour plots found within this report. 


. B. AUTOMATION 

A subroutine, CONISD, was found in the computer program 
library that was designed to produce a contour map of irregu- 
larly spread data points. Each data point is a triad of x, y 


and z values where 


NEQU 


This procedure is nicely suited for plotting the flux contours 


required since the flux per unit area has been calculated at 
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each increment in range and scattering angle. In transforming 
evenly spaced range and angle intervals from polar to cartesian 
coordinates, irregularly spaced intervals are formed. By 
inspection the flux varies quite linearly in the range direction. 


The region to be contoured is subdivided into triangular 


areas, using acute triangles as much as possible. The method 
of triangularization is based on that of Ref. 30. The list of 
triangles is then analyzed for adjacencies. For a given con- 


tour level, each outer boundary of the area to be contoured 

is checked for a possible entry value. Upon finding a value, 
the contour line is traced from triangle to triangle until it 
exits the area. As the line is traced from segment to segment, 
the four nearest values are linearly intérpolated and the 
resulting value stored. After all lines have been found and 
intersections stored, another subroutine is called to. fit 
smooth curves to the stored values. The resulting curves are 


the desired contours of equal photon flux. 


“2 ADDITIONS TO DRLITE 

A section was created within DRLITE to drive the subroutine 
CONISD after computation of relative flux at desired grid 
points. Inputs required by CONISD are the number of grid points, 
each corresponding triad of values, the value of each contour 
desired, scaling factors for plotting the output and designa- 
tion of which data points are to be boundary points. 

Because data used to generate the contours was generated by 


a Monte Carlo routine, a measure of statistical significance had 
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to be included before actual plotting. The measure used was 

based purely on the number of photon crossings at each parti- 
cular angular bin at every range. If the number of crossings 
was less than three, the corresponding data triad was removed 


from the list to be used in generation of contours. 


D. DISCUSSION 

The automation method explained above, although useful in 
reducing manual labor, is quite consumptive of CPU time.  Com- 
puter time goes up almost expotentially with the number of 
grip points used in the routine. This can be attributed to 
the number of searches necessary in passing from segment to 
segment. Although the general idea behind automation of the 
contour plotting is good, the actual implementation of such 
methods must be done with caution. 

As an aside, the subroutine CONISD also has the capability 
to have internal cut out boundaries placed into its logic. 
Future application of this feature may be to draw equal photon 


flux contours for light propagation around defined objects. 
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APPENDIX E 
I. DERIVATION, DOCUMENTATION AND VERIFICATION 
OE EEPEEEUVSDOATIENUATTON COEPBPTOCTENT METHOD 
A. INTRODUCTION 
Gordon [9] presents the concept of an effective attenuation 

coefficient (EAC) which considerably simplifies multiple scat- 
tering calculations. Concise expressions for the total flux 
through an aperture and the beam-spread function are derived 
in terms of the EAC. A program was written to evaluate these 
expressions so that comparison with other light scattering 
models could be accomplished. This Appendix derives the 
expressions given by Gordon, documents the routine written to 
evaluate the expressions and shows verification of program 


accuracy. 


B. DERIVATION 

Gordon first suspected that simple formulas might describe 
multiple scattering after examining Duntley's [29] measurements 
of the fraction of power emitted from a collimated source which 
reaches a circular aperture subtending a cone of half-angle 9 
when viewed from the source and located a distance r from it. 
Gordon noticed that on semi-log paper the data formed straight 
lines to about 15 extinction lengths when plotted as a function 
of range. The expression for flux reaching the aperture where 


Р іѕ the source strength can be expressed as [9]: 


TI 





e 9e (OR 


F(8,R) = Fo (1) 


where a,(8) is the EAC as a function of 9. 

Consider the geometry of Figure 38. A unidirectional point 
source is centrally located on the plane x = O. On the plane 
x = Ris located a circular aperture which subtends a cone of 
half-angle 8 with the point source. The total flux through 


the aperture is given by [9]: 


2n -9' 


1 
Р = F exp 63 TI Po) dat «con (2) 
o `о о 


where, 


do - sin 9 d0 dÜ - differential solid angle 


S - scattering coefficient in km | 
æ = extinction coefficient in kn i 
t = x/R = normalized range 


"T -l[tan Ө 
е (тк 


and P(8) is the normalized phase function, thus it satisfies the 


relation, 


2т oq 
If P($) ао = 1 (3) 
o~ o 
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Figure 38 


Geometry of EAC Method 
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A simple physical description of equation 2 is that it 
accounts for flux absorbed and scattered from the beam in the 
first term of the exponent and also accounts for flux that is 
scattered back through the aperture in the integral term of 
the exponent. 

Comparing equation 1 with equation 2 the effective attenu- 
ation coefficient is found to be 


nenn: 
Q O (9) = (esf fme d (4) 
а 
O 


In adapting this equation to machine computation, a normali- 
zed phase function is necessary. A reasonable choice is the 


Henyey-Greenstein phase function [7], 


Cres) 


роз = eo e s 
4-T(1*G?-2Gcos0)! ? 


-1<6<1 (5) 
Because the Monte Carlo model makes use of this phase function 
also, comparison of the two models is possible. 

Now that the flux through an aperture is known, it can be 
used to derive the irradiance caused by a unit strength uni- 
directional point source at each point of a selected target 
plane. This irradiance distribution is called the beam spread 
Function, H„(8,R). Let the initial point source have unit 
strength. By considering the flux through a small annulus of 


width R d8 located at (R,9), H„(8,R) can be written as [89], 
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L dE 


Hg(9.R) = тзт 48 > | 
where dF/d9 is the derivative of equation 2 with A = 14 To 
evaluate this derivative, Leibnitz rule is used [33]. It states 
that 1f 

S4 CO 
ECT) -| D(x,t)dt (7) 
y (4) 


where a(t) and b(t) are differentiable functions of t and where 


ß(x,t) and IL are continuous in x and t, then 
pt) 


dr | 3D OCA ах+й[ъ(+) ,+1995® ate) 12 (8) 


dt 
Жез, 





Take the derivative of equation 2, 


] cn O 


ПЕ а S 
ag = Fag (feel): (9) 
- JO О 
8! 
2TFSR [| P(@)sin Ө «| dt (10) 
О 


and apply Leibnitz rule where 


8! 


й(+,0) Ji Р(9)ѕіпе аб а(+) = 0 (11) 
О 


HH il 
= 


b(t) 


it is now apparent that the derivative becomes, 


3-16 





! 
dr _ 2TFSR 52 [ec sind dd} dt (12) 
do | 85 


The derivative that remains is of an integral whose range is 
over the same independent variable, 0, therefore the result is 
the product of the integrand evaluated at the upper limit and 


the derivative of the upper limit, i.e., 


oe sinodo s P(O') sim(O') (01) (13) 


but 





l-t 


9 -1 /tan9 І 
= | tan = Se  — (14) 
I | ( =) Foy (1-t)cos?0 


SO 1t is clear that 





1 
c = arrer | P(8')sin(0') ——— CIE (15) 
O Pe J|e me 
Using the small angle approximation, tan9 = 8 = sin0 and letting 
s -if[ 8 
u = tan (rex) (16) 
A — (17) 


Pe" (1-t) 


equation 15 takes the form 
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1/2 
= = ита] P(u) cos(u) du (18) 
O 


Using equation 18 in equation 6 the beam spread function is 


ПЕ? 
H,(9,R) = =s j P(u) cosu du (19) 
O 


RO 

A useful calculation in comparing this beam spread function 
with that found using the Monte Carlo method is taking the 
relative logarithm of the beam spread function. 
©. DADAPTATTON OE EPEECTIVE ATTENTUATION COEFFICIENT METHOD TO 

MACHINE COMPUTATION 

The Effective Attentuation Coefficient method lends itself 
nicely to machine computation. The inputs required by the 
program are the scattering coefficient, the extinction coef- 
`- ficient and the Henyey-Greenstein G parameter. Also entered 
are the initial range and theta, the increment in range and 
theta and the number of grid points desired in each dimension. 

The routine written is called EAC.  EAC calls on a library 
subroutine DQG32 which integrates functions FCC and FCT inside 
of which the integrands of equation 4 and 19 are defined, 
respectively. DQG32 uses 32 point Gaussian quadrature which 
integrates polynomials up to degree 63 exactly. 

The output of EAC can be received by numerical matrix tabu- 
lation or by graphical contour plotting as described in Appendix 


D. Examples of these plots are found in Figures 6 through 9. 
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D. VERIFICATION OF EAC 

The results from EAC, the NPS routine, were compared to 
figures obtained by Gordon [9] using the same method. There is 
a small dependence on phase function for small optical thick- 
nesses for which this method is applicable. Thus in comparing 
results with Gordon, a slight deviation is noticed due to the 
fact that the phase function used may not have been exactly the 
same. A Henyey-Greenstein G parameter of .94 was used to gen- 
erate the figures presented here. Figure 39 compares the 
results of EAC with Gordon for irradiance as a function of Ө 
at a range of 2.12 extinction lengths. Agreement is seen to be 
quite good. Figure 40 compares the results of EAC with Gordon 
for irradiance as a function of range with a scatter angle of 
10 degrees. Agreement is good here also. These are a few of 
many checks on routine accuracy which lead to much confidence 


in the routine. 
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APPENDIX F 


l. CHECKS ON POSSIBLE ERRORS 


A. GENERAL 

Each computational procedure employed in this simulation 
Was verified when possible and every effort was made to ensure 
the entire simulation functioned properly. Results of various 
theories were compared using identical parameters to build 
confidence at intermediate steps along the path toward final 
results. Procedures used to document the correct behavior 
of the computer routines used in this simulation are described 


in this Appendix. 


Be oracle tC CHECKS 

One of the secondary objectives of this report was to com- 
pare results of different theories in predicting light transfer 
through a scattering absorbing medium. The ability to compare 
results of separate theories stands by itself as a verification 
of the accuracy of each theory. 

Reference 6 states numerous initial methods used to verify 
the mechanics of the computer simulation in its first stages. 
In converting the routine from one that simulates an infinite 
homogeneous medium to one that contains a cloud, additional 
checks were required. Many photons were traced manually by 
hand calculations to ensure that all nine possible combinations 


of transfer across boundaries were correctly handled by the 
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routine. ES were found and corrected on many an occasion. 
In the situation where the transmission parameters and the 
phase function of the cloud were the same as those of the 
surrounding medium, one would expect the results to be identi- 
cal to those found when simulating a homogeneous medium. This 
in fact was the case. Results were identical to well within 
Statistical scatter. 

Changeover from measurements in terms of extinction lengths 
to measurements in terms of real dimensions also required brief 
checking. One would expect the results to be identical when 
comparing a medium with an extinction coefficient of one in- 
verse kilometer to a medium described in terms of optical 
thickness. This is so because one unit of distance is the 
Same in both cases, one kilometer. This is the case. In com- 
paring trials where the first has a distance between accounta- 
bility shells of twice the second trial but an extinction 
coefficient of half the second trial with the albedo the same, 
the number of photons crossing each bin should be the same. 
Once again this is the case. 

In some cases the results were checked by intuition only 
as no previous work with similar presentation was found. 
Specifically, the equal photon flux contours when passing 
through clouds were checked only by what one would expect them 
to look like. 

Many test runs were made using the model generating a 


scatter angle weighted by an arbitrary phase function. The 
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angles generated are in fact weighted by areas of panels and 
linearly within, panels. 

Confidence in the statistical nature of a Monte Carlo 
routine is related to the number of photon histories used in 


each trial. The smallest number of photons possible were used 

to create outputs of the desired accuracy in order to ninimise 
computer time consumption. In most cases of spatial spread 
thousands of photon histories were generated. In time spread 
cases as few as 100 were required at times. 


Other checks made are found in Ref. 6 as well as within 


the body of this work. 
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APPENDIX 6 
I. DERIVATION OF CLOSED FORM EXPRESSION 
FOR TIME SPREAD 


Stotts [26] derives a closed form expression for time spread 
as follows. From vector analysis, the average pathlength of a 


photon per unit time is 


GR _ |.d2,2,.dy,2,.42,2|% 
dt ES rr «ib 


which reduces easily, where г = хі = yj, to 


dR dr); 21% 
dE = E £e dz 


From basic scattering theory [31] the projected angle of scatter- 


ing is approximated by 
27% 
[aot Ye 
where wis the albedo of single scattering, 1 is the optical 
thickness of the scattering region and y, 1s the RMS scatter 
angle. Thus the projected mean-squared transverse displacement 


is 


= 1 
lär| = z(wot v2) 
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thus 


and 


7 
ПЕ = (142 wT Sé dz 


Integrating results in 


_ 8z 3 сурт 
т 270 T =| (5 с B f 


The average multipath time spread is defined as the difference 
between the average transient time incurred from multiple 
scattering and the normal transient time in the absence of 


scattering. Thus, Stotts has for time spread, 


2 ‚30 9 2 115 
= + = -—À 
E = zi os б т ү J | 18 
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